
// Autor: Andre  Ricardo Goncalves
// Dia: 03/04/2008



// Arquivo que contem a implementacao dos Metodos Diretos e Indiretos

// METODOS DIRETOS:
//                  Eliminacao de Gauss c/ Pivoteamento Parcial
//                  Decomposicao LU
//                  Cholesky

//  METODOS INDIRETOS
//                  Jacobi
//                  Gauss-Seidel



// ***************************************************************************************
//              FUNCOES PARA O METODO DE ELIMINACAO DE GAUSS C/ PIVOTEAMENTO
// ***************************************************************************************

// funcao que chama as funcoes necessarias no metodo de Gauss
function [vet_solucao] = main_gauss_pivot( a , b )

// verifica se o determinante da matriz de coeficientes e' igual a zero
if ( det(a) ==  0 ),
    // caso for, este sistema nao pode ser resolvido por este metodo
    vet_solucao = 'Determinante igual a 0. Infinitas Solucoes. Nao pode ser utilizado o metodo de Gauss c/ pivotamento.';
    break;
    
else
  printf("Matriz de coeficientes:");
  disp(a);
  printf("\nTermos independentes:");
  disp(b);

  // chama a funcao para realizar o processo de triangularizacao
  mat_triang = triangulariza( a , b );

  if( type(mat_triang) == 10 ),
       vet_solucao = mat_triang; 
       break;
   else
  	// chama a funcao para realizar o processo de retrosubstituicao
  	vet_solucao = retrossub( mat_triang );
  end
 
  // retorna o vetor solucao
end

endfunction


// Funcao que realiza a tringularizacao
// da matriz expandida ( matrizCoef | b )
function  [mat_triang] = triangulariza( matrizCoef , b )

// verifica se o determinante da matriz de coeficientes e' negativa
// caso for o sistema nao tera solucao

// cria a matriz expandida
matriz_exp = [ matrizCoef b ];

//mat_triang = null;

// Varre todas as linhas da matriz expandida
// realizando o processo de triangularizacao
for i=1:size(matriz_exp,'r')-1,

     // Caso o encontre algum pivo igual a 0 (zero), deve buscar uma linha
     // abaixo que nao contem zero, isto evita a tentiva de divisao por 0 (zero)
     // numero de linhas da matriz expandida
     nlinhas = size(matriz_exp , 'r');

     // obtem o maior valor absoluto e a posicao do elemento no vetor linha abaixo da linha pivo 
     [ m , n ] = maxi( abs( matriz_exp(i:nlinhas,i) ) );

     // nao precisa trocar, caso o proprio elemento pivo seja o maior 
     if( n(1) <> 1 ),
     		aux = matriz_exp( i , : );
     		matriz_exp( i , : ) = matriz_exp( n(1)+(i-1) , : );
     		matriz_exp( n(1)+(i-1) , : ) = aux;
     end

     if( m == 0 ),
         printf("Matriz Pivoteada: ");
         disp(matriz_exp);
         mat_triang = "Nao existe mais nenhum pivo diferente de 0. Este metodo nao eh adequado, utilize outro metodo.";
         return mat_triang;
     end

     // calcula os novos valores que serao atribuidos as linhas inferiores a linha pivo
     for j=i+1:size(matriz_exp,'r'),
           m = matriz_exp(j,i) / matriz_exp(i,i);
           matriz_exp(j,:) =  matriz_exp(j,:) - matriz_exp(i,:) * m ;
     end;
end

// retorna a matriz triangularizada
mat_triang = matriz_exp;

endfunction


// Funcao que realiza a substituicao dos valores das variaveis
// Processo chamado de retrosubstituicao
function [vetor_solucao] = retrossub( A , b )

  mat_exp = [ A b ];
  nlin = size(mat_exp,'r');
  
  // atribui para o vetor de solucao o unico valor conhecido ate o momento
  // que e' o xn
  xn = mat_exp(nlin,nlin+1)/mat_exp(nlin,nlin);
  L = nlin-1;
  
  // cria um vetor solucao, para ir armazenando os resultados obtidos pela
  // retrosubstituicao 
  sol = zeros(nlin,1);  sol( nlin ) = xn;

  // Varre de baixo para cima todas linhas, substuindo os valores de Xn ja conhecidos 
  // e calculando os Xn-1... X1, ainda nao conhecidos.
  while L >= 1,
       sol(L) = ( mat_exp(L,nlin+1) - sum(mat_exp(L,L+1:nlin) .* sol(L+1:nlin)' ) ) / mat_exp(L,L);
       L = L - 1;
  end;

  // retorna o vetor solucao
  vetor_solucao = sol;

endfunction







// ***************************************************************************************
//              FUNCOES PARA O METODO DA DECOMPOSICAO LU
// ***************************************************************************************


// funcao que chama as funcoes necessarias no metodo de Gauss
function [vet_solucao] = main_lu( a , b )

// verifica se o determinante da matriz de coeficientes e' igual a zero
if ( det(a) ==  0 ),
    // caso for, este sistema nao pode ser resolvido por este metodo
    vet_solucao = 'Determinante igual a 0. Infinitas Solucoes. Nao pode ser utilizado o metodo da Decomp. LU.';
    break;

// verifica  se a matriz de coeficientes e' uma matriz definida positivamente
elseif( def_positivo(a) == 1 )
    vet_solucao = 'Matriz nao eh definida positivamente.';
    break;

else
  printf("Matriz de coeficientes:");
  disp(a);
  printf("\nTermos independentes:");
  disp(b);

  // chama a funcao para realizar o processo de triangularizacao
  [ L , U ] = Dec_LU( a , b );

  // imprimi as matrizes L e U na tela 
  printf("L:"); disp(L);
  printf("U:"); disp(U);

  // chama a funcao de substituicao para frente para resolver
  // o sistema Ly=b
  y = subst_frente( L , b );

  // chama a funcao para realizar a retrosubstituicao para o sistema Uy=x
  vet_solucao = retrossub( U , y );

  // retorna o vetor solucao
  
end

endfunction

// funcao que realiza a decomposicao matriz da A em LU
function  [ L , U ] = Dec_LU( A , b )

      n = size( A , 'r' );
      U = zeros( n , n );  L = zeros( n , n );

      U(1,:) = A(1,:);

      for i=1:n,
           for j=1:n,
                if( i > j ),
                      L(i,j) = ( A(i,j) - L(i,1:j-1) * U(1:j-1,j) ) / U(j,j);
                else
                      U(i,j) = A(i,j) - L(i,1:i-1) * U(1:i-1,j);
                end
           end
           L( i , i ) = 1;
      end


endfunction




// Funcao que realiza a substituicao dos valores das variaveis
// Processo chamado de substituicao para frente
 function [vetor_solucao] = subst_frente( A , b )

  nlin = size(A,'r');
  // atribui para o vetor de solucao o unico valor conhecido ate o momento
  // que e' o x1 ( primeiro elemento do vetor solucao )
  x1 = b(1,1)/A(1,1);

  // cria um vetor solucao, para ir armazenando os resultados obtidos pela
  // substituicao para frente 
  sol = zeros(nlin,1); 
  sol( 1 ) = x1;

  // Varre de cima para baixo todas linhas, substuindo os valores de X1 ja conhecidos 
  // e calculando os X1, ..., Xn, ainda nao conhecidos.
  for i=2:nlin,
       sol(i) = ( b(i) - sum(A(i,1:nlin) .* sol') )/A(i,i);
  end;

  // retorna o vetor solucao
  vetor_solucao = sol;

endfunction


// funcao que verifica se a matriz de coeficientes
// e' definida positivamente
// Se sim: retorna 0
// Se nao: retorna 1
function [ resp ] = def_positivo( A )

      resp = 0;
      for i=1:size(A,'r'),
          determinante = det( A(1:i,1:i) );
          if( abs(determinante) <= 10^(-6) ),
              resp = 1;
          end;
      end

endfunction






// ***************************************************************************************
//              FUNCOES PARA O METODO DE CHOLESKY
// ***************************************************************************************


// funcao que chama as funcoes necessarias no metodo de Cholesky
function [vet_solucao] = main_cholesky( a , b )

// verifica se o determinante da matriz de coeficientes e' igual a zero
if ( det(a) ==  0 ),
    // caso for, este sistema nao pode ser resolvido por este metodo
    vet_solucao = 'Determinante igual a 0. Infinitas Solucoes. Nao pode ser utilizado o metodo de Cholesky.';
    break;
    
// verifica  se a matriz de coeficientes e' uma matriz definida positivamente
elseif( isSimetrica( a ) == 1 )
    vet_solucao = 'Matriz nao eh simetrica. Utilize outro metodo.';
    break;

elseif( def_positivo_chol(a) == 0 )
    vet_solucao = 'Matriz nao eh definida positivamente.';
    break;

else
  printf("Matriz de coeficientes:");
  disp(a);
  printf("\nTermos independentes:");
  disp(b);

  // chama a funcao para realizar o processo de triangularizacao
  [ U ] = MyCholesky( a , b );

  // imprimi a matriz U na tela 
  printf("Matriz U:"); disp(U);

  // chama a funcao de substituicao para frente para resolver
  // o sistema U'y=b
  y = subst_frente( U' , b );

  // chama a funcao para realizar a retrosubstituicao para o sistema Uy=x
  vet_solucao = retrossub( U , y );

  // retorna o vetor solucao

end

endfunction


// funcao que realiza a decomposicao matriz da A em U
function  [ U ] = MyCholesky( A , b )

      n = size( A , 'r' );
      U = zeros( n , n );
      U(1,1) = sqrt( A(1,1) );
      for i=1:n,
           for j=1:n,
                if( i == j ),
                      U(i,j) = sqrt ( ( A(i,j) - sum( (U(1:i-1,i))^2 ) ) );
                elseif( j > i)
  		                 U(i,j) = ( A(i,j) - ( U(1:i-1,i) * U(1:i-1,j) ) ) / U(i,i);
                end
           end
      end

endfunction


// Funcao que verifica se uma matriz eh simetrica
// Se simetrica : 0
// Se NAO-simetrica : 1 
function [ resp ] = isSimetrica( A )

     if( isequal( A , A') == %F ), // nao eh simetrica
          resp = 1;
     else                         // eh simetrica
          resp = 0;
     end

endfunction


// funcao que verifica se a matriz de coeficientes
// e' definida positivamente
// Se sim: retorna 1
// Se nao: retorna 0
function [ resp ] = def_positivo_chol( A )

      resp = 1;
      for i=1:size(A,'r'),
          determinante = det( A(1:i,1:i) );
          if( determinante < 10^(-6) ),
              resp = 0;
          end;
      end

endfunction




// ***************************************************************************************
//              FUNCOES PARA O METODO DE JACOBI
// ***************************************************************************************


function [ vetorSolucao ] =  main_jacobi( A , b )

if( det(A) == 0 ),
   vetorSolucao = "Determinante igual a 0. Infinitas Solucoes. Nao pode ser utilizado o metodo de Jacobi.";
   break;
else
   erro = 1*10^(-6);   // erro aceitavel
   niteracoes = 250;
   
   printf("Matriz de coeficientes:");
   disp(A);
   printf("\nTermos independentes:");
   disp(b);
   printf("\n\nErro Aceitavel: %f \n",erro);
   printf("Numero Maximo de Iteracoes: %d \n",niteracoes);
   
   // verifica se existe elementos nulos na diagonal principal
   // caso existir tenta retira-lo utilizando o pivoteamento parcial
   m_exp = [ A b ];
   m_exp = pivoteamento( A , b );

   if( type( m_exp ) == 10 ),
        vetorSolucao = m_exp;
        return vetorSolucao;
   end

   A = m_exp( : , 1:size(m_exp,'r'));
   b = m_exp( : , size(m_exp,'r')+1);

   if( estr_diag_dom(A) == 1),
      disp('O sistema nao eh estritamente diagonalmente dominante. O metodo de Jacobi nao garante a convergencia.');
      halt();
   end

   [ vetorSolucao , iteracoes ] = jacobi( A , b , erro , niteracoes);

end

endfunction




// Implementa o metodo de Jacobi
// A: matriz de coeficientes
// b : termos independentes
function [ x , it ] = jacobi( A , b , erro_aceit , niteracoes )

       sol = b; // solucao inicial
       erro_rel = 1;
       it = 0;

       while( (erro_aceit < erro_rel) & it < niteracoes ),

           for i=1:size(A,'r'),
                x(i) = (-1)*( (A(i,:)*sol-( A(i,i)*sol(i)+b(i) ) )/ A(i,i));
           end
           erro_rel = abs((norm(sol)-norm(x))) / norm(sol);
           sol = x;
           it = it + 1;
	   if( it == niteracoes ),
	       x = string("O metodo nao convergiu para ")+string(it)+string(" iteracoes.");
	       return;
	   end 

       end
       x = sol;
       
endfunction


// Funcao que realiza o processo de pivoteamento
// utilizada nos metodos iterativos do jacobi e do gauss-seidel
function [ mat_triang ] = pivoteamento( A, b )


    matriz_exp = [ A b ]; 
    nlinhas = size(matriz_exp , 'r');
    
    for i=1:size(matriz_exp,'r')-1,
      
         inc_linha =  1;

          if( matriz_exp(i , i) == 0 ),
                     // Caso o encontre algum pivo igual a 0 (zero), deve buscar uma linha
                     // abaixo que nao contem zero, isto evita a tentiva de divisao por 0 (zero)
                     // numero de linhas da matriz expandida


                     // obtem o maior valor absoluto e a posicao do elemento no vetor linha abaixo da linha pivo 
                     [ m , n ] = maxi( abs( matriz_exp(i:nlinhas,i) ) );

                     if( m == 0 ),
                         printf("Matriz Pivoteada: ");
                         disp(matriz_exp);
                         mat_triang = "Nao foi possivel retirar o elemento nulo da diagonal principal. Utilize outro metodo.";
                         return mat_triang;
                     end 
               
                     // nao precisa trocar, caso o proprio elemento pivo seja o maior 
                     if( n(1) <> 1 ),
                     		 aux = matriz_exp( i , : );
                     		 matriz_exp( i , : ) = matriz_exp( n(1)+(i-1) , : );
                     		 matriz_exp( n(1)+(i-1) , : ) = aux;
                     end

           end
   end
  
  if( matriz_exp(  nlinhas , nlinhas ) == 0 ),
           printf("Matriz Pivoteada: ");
           disp(matriz_exp);
           mat_triang = "Nao foi possivel retirar o elemento nulo da diagonal principal. Utilize outro metodo.";
           return mat_triang;
  end
   
  mat_triang = matriz_exp;
  
endfunction



//Funcao que verifica se a matriz eh estritamente diagonalmente dominante
function [ res ] = estr_diag_dom( A )

       res = 0;
       for i=1:size(A,'r'),
		if( A(i,i) < (sum( abs(A(i,:)) )-abs(A(i,i))) ),
		   res = 1;
		end
       end

endfunction


// ***************************************************************************************
//              FUNCOES PARA O METODO DE GAUSS-SEIDEL
// ***************************************************************************************



// Funcao principal do gauss-seidel
// realiza a verificacao dos criterios que o sistema
// deve suportar para que o mesmo possa ser solucionado pelo metodo
// de gauss-seidel
function [ vetorSolucao ] =  main_gauss_seidel( A , b )

// Verifica se o determinante da matriz de coeficientes eh nulos
// caso for, o sistema nao pode ser resolvido pelo metodo de gauss-seidel
if( det(A) == 0 )
   vetorSolucao = "Determinante igual a 0. Infinitas Solucoes. Nao pode ser utilizado o metodo de Gauss-Seidel.";
   break;
else
   // definindo o erro aceitavel e o numero maximo de iteracoes
   // que sera suportado
   erro = 1*10^(-6);   // erro aceitavel
   niteracoes = 250;
   
   printf("Matriz de coeficientes:");
   disp(A);
   printf("\nTermos independentes:");
   disp(b);
   printf("\n\nErro Aceitavel: %f \n",erro);
   printf("Numero Maximo de Iteracoes: %d \n",niteracoes);

   // verifica se existe elementos nulos na diagonal principal
   // caso existir tenta retira-lo utilizando o pivoteamento parcial
   m_exp = [ A b ];
   m_exp = pivoteamento( A , b );

   // verifica que o pivoteamento nao retornou algum erro
   if( type( m_exp ) == 10 ),
        vetorSolucao = m_exp;
        return vetorSolucao;
   end

   A = m_exp( : , 1:size(m_exp,'r'));
   b = m_exp( : , size(m_exp,'r')+1);

   // verificando a convergencia pelo criterio de sassenfeld
   if( criterio_sassenfeld(A) == 1 ),
    disp("O sistema nao satisfaz o criterio de sassenfeld. O metodo de Gauss-Seidel nao garante a convergencia.");
      halt();
   end
   // verificando a convergencia pelo criterio da matriz de coeficientes ser estritamente diagonalmente dominante
   if( estr_diag_dom(A) == 1),
   disp("O sistema nao eh estritamente diagonalmente dominante. O metodo de Gauss-Seidel nao garante a convergencia.");
      halt();
   end
   // chama a funcao que realiza o metodo de gauss-seidel  
   [ vetorSolucao , iteracoes ] = gauss_seidel( A , b , erro , niteracoes);

end

endfunction




// Implementa o metodo de Gauss-Seidel
// A: matriz de coeficientes
// b : termos independentes
function [ x , it ] = gauss_seidel( A , b , erro_aceit , niteracoes )

       sol = b; // solucao inicial
       erro_rel = 1;
       it = 0;

       // enquanto o criterio de parada nao for satisfeito
       // 1. Erro aceitavel for menor do que o erro de aproximacao
       // 2. Um numero de iteracoes limite nao foi satisfeito
       while( (erro_aceit < erro_rel) & it < niteracoes ),

	   // guarda a solucao anterior para calcular o erro relativo
           anterior = sol;
           // metodo de gauss-seidel
           for i=1:size(A,'r'),
                sol(i) = (-1)*( (A(i,:)*sol-( A(i,i)*sol(i)+b(i) ) )/ A(i,i));
           end
           // calculo do erro relativo
           // para determinar um valor escalar de um vetor 
           // utilizou-se o calculo da norma
           erro_rel = abs((norm(anterior)-norm(sol))) / norm(anterior);
           it = it + 1;

           // se o numero de iteracoes foi satisfeito 
           // e o erro nao foi menor do que o erro aceito, entao
           // termina o processo e afirma que o metodo nao convergiu
	   if( it == niteracoes ),
	       x = string("O metodo nao convergiu para ")+string(it)+string(" iteracoes.");
	       return;
	   end 

       end
       x = sol;
endfunction


// Funcao que implementa o criterio de sassenfeld
// Retorna 0 se satisfaz
// Retorna 1 se NAO satisfaz
function [ res ] = criterio_sassenfeld( A )

         // inicilizo o vetor Beta com Zeros 
         B = zeros( 1 , size(A,'r') );
         n = size(A,'r');

	      // criando o veto Beta 
         for i=1:size(A,'r'),
                a_estrela = A(i,:)/A(i,i);
    	 	         B(i) =  abs(a_estrela(1:i))*B(1:i)' + sum(abs(a_estrela(i+1:n)));
         end

	       // calcula o maximo elemento de beta
         maxBeta  = maxi(B);

         // satisfaz o criterio de sassenfeld
         if( maxBeta < 1 ),
             res = 0;
         else  // nao satisfaz o criterio de sassenfeld
             res = 1;
         end

endfunction







// ***********************************************************************************************
//             FUNCAO QUE CHAMA OS METODOS PARA SEREM EXECUTADOS COM O MESMO TESTE
// ***********************************************************************************************
  
// Metodos de  Resolucao de Sistemas Lineares
  
function  chama_metodos()  
           
            // abre a janela pra escolha da matriz com os coeficientes e termos independentes.
            nome_arq = tk_getfile("*.txt",Title="Clique no arquivo com a matriz expandida.");
            
            // carrega a matriz com mgetl
            // o mgetl apens enxegar a matriz como sendo linhas de strings, ele nao
            // consegue distinguir colunas, portanto utilizamos o mgetl para obter o numero de linha da matriz
            // como estamos carregando a matriz expandida, sabemos que sua dimensao e' n x n+1, sabendo o numero
            // de linhas, pelo mgetl, consequentemente saberemos o numero de colunas.
            // assim podemos utilizar o read(), que consegue enxer uma matriz, com linhas e colunas
            mat = mgetl(nome_arq,-1);
            
            // calculando o tamanho da matriz_expandida
            numLinhas = size(mat,'r');
            ncol = numLinhas + 1;
            
            // obtendo a matriz do arquivo com o read()
            matriz = read( nome_arq , -1 , ncol ); 
            
            // limpar a tela
            clc;
            
            // chamando a funcao que executa o metodo de gauss
            printf("--------------------------------------------------------\n")
            printf("::    METODO ELIMINACAO DE GAUSS C/ PIVOTEAMENTO      ::\n");
            printf("--------------------------------------------------------\n\n")

            vet = main_gauss_pivot( matriz( : , 1:ncol-1 ) , matriz( : , ncol ) );
            printf("Resposta do Metodo: \n");
            disp(vet);
            halt();
            clc;  
            
            // chamando a funcao que executa o metodo da Decomposicao LU
            printf("--------------------------------------------------------\n")
            printf("::           METODO DA DECOMPOSICAO LU                ::\n");  
            printf("--------------------------------------------------------\n\n")  
             
            vet = main_lu( matriz( : , 1:ncol-1 ) , matriz( : , ncol ) );
            printf("Resposta do Metodo: \n");
            disp(vet);
            halt();
            clc;
            
           // chamando a funcao que executa o metodo de Cholesky
            printf("---------------------------------------------------------\n")
            printf("::                 METODO DE CHOLESKY                  ::\n");
            printf("---------------------------------------------------------\n\n")
            vet = main_cholesky( matriz( : , 1:ncol-1 ) , matriz( : , ncol ) );
            printf("Resposta do Metodo: \n");
            disp(vet);
            halt();
            clc;
            
           // chamando a funcao que executa o metodo de Jacobi
            printf("--------------------------------------------------------\n")
            printf("::                  METODO DE JACOBI                  ::\n");
            printf("--------------------------------------------------------\n\n")
           
            vet = main_jacobi( matriz( : , 1:ncol-1 ) , matriz( : , ncol ) );
            printf("Resposta do Metodo: \n");            
            disp(vet);
            halt();
            clc;
            
             // chamando a funcao que executa o metodo de Gauss-Seidel
            printf("--------------------------------------------------------\n")
            printf("::             METODO DE GAUSS-SEIDEL                 ::\n");
            printf("--------------------------------------------------------\n\n")
            
            vet = main_gauss_seidel( matriz( : , 1:ncol-1 ) , matriz( : , ncol ) );
            printf("Resposta do Metodo: \n");            
            disp(vet);
            halt();
            clc;
            disp('Fim da execucao dos metodos');            
  

endfunction
